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RELATIVE ADVANTAGES OF THIN-LAYER 
NAVIER-STOKES AND INTERACTIVE 
BOUNDARY LAYER PROCEDURES* 

SUMMARY 


Numerical procedures for solving the thin-shear-layer Navier-Stokes equations and for the 
interaction of solutions to inviscid and boundary-layer equations are described and evaluated. To 
allow appraisal of the numerical and fluid dynamic abilities of the two schemes, they have been 
applied to one airfoil as a function of angle of attack at two slightly different Reynolds numbers. 
The NACA 0012 airfoil has been chosen because it allows comparison with measured lift, drag, 
and moment and with surface-pressure distributions. Calculations have been performed with 
algebraic eddy-viscosity formulations, and they include consideration of transition. The results 
are presented in a form that allows easy appraisal of the accuracy of both procedures and of 
the relative costs. The interactive procedure is computationally efficient but restrictive relative 
to the thin-layer Navier-Stokes procedure. The latter procedure does a better job of predicting 
drag than does the former. In both procedures, the location of transition is crucial for accurate 
or detailed computations, particularly at high angles of attack. When the upstream influence of 
pressure field through the shear layer is important, the thin-layer Navier-Stokes procedure has an 
edge over the interactive procedure. 


INTRODUCTION 


It is generally accepted that the Navier-Stokes equations correctly represent fluid-flow phe- 
nomena. Since the unsteady, three-dimensional equations can usually be solved for flows in which 
small-scale fluctuations are unimportant, emphasis has been placed on particular reduced forms 
such as those appropriate to regions of inviscid flow and boundary layers. In recent years, and 
with the application of numerical solution procedures in mind, attention has also been paid to 
the Reynolds-averaged Navier-Stokes equations and to various further-reduced forms, including 
their so-called parabolized forms and the thin-layer Navier-Stokes (TLNS) equations. 

In this report, we are concerned with calculating the properties of steady incompressible flows 
around airfoils at angles of attack from 0° up to values where separation can occur near the trail- 
ing edge of the suction surface. The consideration of airfoils itself introduces a simplification in 
that the corresponding Reynolds-averaged equations have only two independent variables. Three 
main approaches to problems of this type are described in the literature and are characterized by 
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different forms of equations. The first possibility is to reduce the Navier-Stokes equations to a 
Reynolds-averaged, two-dimensional form resulting in elliptic equations that can be numerically 
solved. The second possibility is to further reduce the equations by neglecting the longitudinal 
diffusion terms leading to the TLNS equations, which are elliptic-parabolic and strongly coupled. 
The third possibility is to solve inviscid-flow equations and to use the resulting pressure distri- 
bution as a condition for the solution of the boundary-layer equations, which are parabolic, with 
or without a second momentum equation to represent any normal pressure gradient. This classic 
approach has been improved in recent years by interaction of the two solutions to ensure that the 
viscous flow is allowed to influence the inviscid flow and vice versa. 

The design requirement that calculations be performed for several geometries and angles of 
attack and for the necessary extension to three dimensions, together with the present emphasis on 
flows with small or no separation, suggests that the choice of methods lies between those methods 
associated with interactive boundary layers (IBL) and those associated with TLNS equations. 
As a consequence, we have performed calculations with the IBL procedure of Cebeci et al. (ref. 
1) and with the TLNS procedure developed at Ames Research Center (refs. 2-4). The results 
to date correspond to a NACA airfoil at several angles of attack. It is expected that these and 
similar results for a range of airfoils will permit the determination of a range of parameters for 
which one or the other method is preferable. 

In this report, we present results that show the influence of transition location and of the 
turbulence model, in particular, models based on the Cebeci-Smith (ref. 5) and Baldwin-Lomax 
(ref. 2) eddy-viscosity formulations. The former formulation has been extensively used in previous 
boundary-layer calculations and the latter in previous applications of the TLNS equations. Since, 
in principle, the Baldwin-Lomax model does not apply to the far-wake region, the model developed 
by Chang et al. (ref. 6) has also been used for wake calculations. 

The TLNS procedure and the IBL procedure are described in the following two sections, which 
are followed by a description of the turbulence models and transition criteria. Computed loads 
acting on a NACA 0012 airfoil are then compared with each other and with experimental data. 
In addition, some details of the flow field, in the form of pressure distributions and displacement- 
thickness distributions, are presented. The report ends with a discussion of upstream influences 
and concluding remarks. 

The authors gratefully acknowledge the assistance provided by Dr. Thomas Pulliam in the 
use of the ARC2D code. 


THIN-LAYER NAVIER-STOKES METHOD 


Some viscous-flow problems can be represented by a set of equations that falls between the 
boundary-layer equations and the full Navier-Stokes equations. This intermediate set of equations 
is applicable to both inviscid and viscous regions, in common with the Navier-Stokes equations; 
therefore, it can be used to compute strong interactions between the two regions. Since they 
couple the viscous and inviscid approximations to a fluid flow, they have been referred to as 
composite equations (ref. 7). A main feature of these equations is the presence of a nonzero, 


2 



normal-pressure gradient, which is necessary to couple and solve simultaneously the viscous and 
inviscid regions. Relative to the Navier-Stokes equations, these composite equations require less 
computational effort because they contain fewer terms. 

The composite, thin-layer equations are generally referred to in the literature as TLNS equa- 
tions. They are obtained by neglecting all streamwise and spanwise derivatives of the viscous 
and turbulence stress, conductive heat-flux terms, and any term involving mixed derivatives. 
These approximations are justified either by physical order of magnitude analysis or by a com- 
putational accuracy argument which amounts- to the observation that since the neglected terms 
cannot be computed correctly with the available grid resolution anyway, why keep them? The 
TLNS equations are primarily based on this latter argument, which is computational rather than 
physical. 

The form of the TLNS equations generally used elsewhere and here, does not satisfy relation- 
ships between metric coefficients in diffusive and conduction terms (ref. 8). The error introduced 
by this oversight is usually insignificant, except when the effective viscosity is relatively large. 
Another limitation of these equations is that the longitudinal-curvature diffusive terms are ne- 
glected. This limitation is a consequence of using the Cartesian velocity components instead of 
curvilinear velocity components with curvilinear coordinates. This problem has been discussed 
in detail by Blottner (ref. 9). 

An implicit numerical method is used to solve the TLNS equations. Some parts of this method 
are essentially those given by Beam and Warming (ref. 3). Other parts are from the numerical 
schemes presented by Steger (ref. 4), Pulliam and Chaussee (ref. 10), Salas et al. (ref. 11), and 
Jameson et al. (ref. 12) Since only steady-state computations are of interest, a diagonal form for 
the Euler equations and a spacially varying time-step are used. These different parts have been 
incorporated in the ARC2D computer code, which is a continuously evolving code. There is no 
detailed documentation of this code. The version 142 of ARC2D code (ARC2D142) was utilized 
with a few crucial modifications for the present study. (The parts of ARC2D142 that have been 
used here are identical to those in ARC2D150, which is the current version of this code.) The 
turbulence modeling modifications are identified in section 4, and the numerical modifications are 
incorporated in the numerical method presented in the appendix. This method is outlined for 
completeness and for documenting a novel derivation of the basic numerical scheme. 


INTERACTIVE BOUNDARY-LAYER METHOD 


The IBL method of Cebeci et al. (ref. 1) has as its essential components Halsey’s inviscid 
procedure (ref. 13), based on the conformal mapping and Fourier analysis techniques, and an 
inverse finite-difference boundary-layer procedure due to Cebeci et al. (ref. 14) This method 
employs a viscous/inviscid iteration procedure in which viscous calculations are performed over 
the airfoil and wake. After each viscous sweep, the external inviscid solution is recomputed, 
and the whole cycle is repeated until a converged solution is obtained. The inverse boundary- 
layer code incorporates the eddy-viscosity formulation, and is able to compute flows with large 
regions of separation without numerical problems. In regions of reverse flow, it uses the FLARE 
approximation, in which the streamwise convective term is set equal to zero in the recirculating 
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region. 

The boundary-layer calculation is initiated on each surface of the airfoil by a direct solution 
in which the boundary condition at the outer edge follows from matching the viscous velocity 
with the latest computed inviscid surface-velocity distribution. As the calculation proceeds down- 
stream, the boundary-layer algorithm is- switched into an inverse mode, in which the viscous-edge 
velocity is computed as part of the boundary-layer solution. This calculation is accomplished 
by applying an interaction procedure suggested by Veldman (ref. 15) and further developed by 
Cebeci et al. (ref. 1) and Cebeci et al. (ref. 14) An overrelaxation scheme is employed, and 
studies have shown that the rate of convergence is considerably improved when it is combined 
with the interactive viscous calculation. In addition, the solution for each successive angle em- 
ploys the previous solution as the starting solution when calculations are required for a range of 
angles of attack. In this way, the converged solution for each angle of attack requires fewer than 
10 iterations. 

The method for calculating the wake region has some novel features which allow results to 
be obtained at high angles of attack. A sudden jump is introduced at the trailing edge by the 
removal of the no-slip boundary condition on the airfoil surface. To account for this jump, a 
small step, with the size related to chord Reynolds number, is employed in the immediate vicinity 
of the trailing edge. For wake calculations involving reverse-flow regions, an additional iterative 
scheme is employed with the FLARE approximation and is based on the homotopy continuation 
method. Studies by Cebeci et al. (ref. 1) have shown that without this added feature, the 
boundary-layer calculations would break down when significant trailing-edge separation region 
is present. Thus, the trailing-edge velocity profile is modified to correspond to an attached flow 
profile; it thereby allows boundary-layer calculations to be performed at the next downstream 
station. The upstream profiles are gradually modified and the downstream profile is recomputed 
until a solution is obtained for the original separated velocity profile. This computational scheme 
is employed at each wake station for which there is flow reversal. Further details are provided in 
reference 1. 


TUBULENCE MODELS AND TRANSITION CRITERIA 


Two algebraic, eddy-viscosity, turbulence models have been used in the present calculations 
The first, (nt) a , is described by Cebeci et al. (ref. 1); the second, (^t)b, is based on the Baldwin- 
Lomax model (ref. 2). Both models stem from the earlier work of Cebeci and Smith (ref. 5). 
The particular formulations used here are described below. 


Eddy-Viscosity Model ( nt) a 

I 

A zero-equation turbulence model, used in the Reynolds-averaged formulation for wall-bounded 
flows that are steady, is due basically to Cebeci and Smith (ref. 5) and is discussed in some detail 
by Cebeci and Smith (ref. 16). The present version of this model is described by Cebeci et al. 
(ref. 1). This model is briefly outlined next. The model is extended into the wake region of the 
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airfoil as done by Chang et al. (ref. 6). This combination of the eddy-viscosity formulation in 
the boundary layer and the wake is labeled as model [nt) a . The nondimensional eddy viscosity 
is represented by 

Htbi = Repavt') tr (-i k ) a (l) 

In the inner layer, 0 < y < y c , the expressions for v and £ are 


(u) 


inner 



{Winner 


OAy 


1 



(2) 


where, with w denoting conditions at the wall and max designating the maximum value in the 
boundary layer, 


y + = y 


Rep t 


du 

w Qy max 


1 1/2 


The use of the maximum value of |3u/3y|, instead of the wall value, accounts for the effect of the 
pressure gradient. In the outer region, y > y c , the expressions for v and £, with the subscript e 
denoting the edge velocity and 6* the usual definition of displacement thickness, are 


[ v )outer — ttei (^)outer — ^ 


( 3 ) 


The distance y c is defined by the continuity of the two eddy viscosities. 

The value of a in the inner layer is 1.0, and its value in the outer layer, suggested by Simpson 
et al. (ref. 17), is 

0.0168 

« = -srr ( 4 ) 


The parameter F is determined from 


F= <1-0 


(du/dx) 

|d«/dy| J 


with the condition that 


The value of /? depends on 


f 1.9919, * 

\ 0.5867, i 


f F > 1.9919 
if F < 0.5867 


Rt = 


(ndu/dy) x 


{fitbldu/dy) max 

and it is determined from the experimental data of Nakayama (ref. 18), 


( 5 ) 


(6) 


( (1 + Rt)/Rt, l f Ft > 1-0 

( 6/[l + 2Rt{2 - Rt)} otherwise 


If R t < 0, then Rt is taken to be zero. 


( 7 ) 
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In equation (l), 7 * r is a transition factor accounting for the transition region that exists 
between a laminar and turbulent flow. It is determined by 


7tr = 1 - exp 


‘ . f x dx 

-G[X - X tr ) / 

J xtr u e - 


( 8 ) 


where the empirical factor G is 


G = 


Re 2 u z e 

g tr 

1200 


(Reu Ctr x t r ) 


-1 34 


and where the subscript tr designates values at the beginning of transition. 

Further, in equation (1), (lk) a is the Klebanoff intermittency factor, which acounts for inter- 
mittent (laminar-turbulent) nature of the outer part of the boundary layer; it is given by 


(7*) a 


1 + 5.5 


(I)' 


-i 


(9) 


where 6 is the boundary-layer thickness. 

The airfoil wake is divided into two regions, the first of which is near the trailing edge and the 
second of which is farther away. In the near-wake, the flow is adjusting to the sudden elimination 
of the wall-boundary condition at the trailing edge, while asymptotically approaching the far- 
wake condition. Therefore, one would expect that the eddy viscosity would be close to that for 
the boundary layers at the trailing edge, and that it would asymptotically approach the far-wake 
expression. Based on the study dealing with wakes subjected to adverse pressure gradients by 
Chang et al. (ref. 6), this is accomplished with the expression 


Htw = IHf + {ute - lHf)exp 


X ^te 
206* e 


) 


( 10 ) 


where 

Ptf — 0.064p7fc {maximum , (Hi 

Here the suscripts l and u denote the lower wake and the upper wake, respectively; and the 
displacement thickness is measured from the location of minimum velocity. 


Eddy- Viscosity Model (//*)*, 

The second zero-equation model is patterned after that of Baldwin and Lomax (ref. 2). The 
Baldwin-Lomax model is patterned after that of Cebeci and Smith (ref. 5). The Cebeci-Smith 
model is difficult to use with Navier-Stokes equations because of the necessity for determining 
the displacement thickness. The Baldwin-Lomax model does not require the location of the outer 
edge of a thin-shear layer. It uses the distribution of vorticity to determine the length scale in 
the outer region of the shear layer. Consequently, it also uses vorticity in the inner region. 


6 



The nondimensional eddy viscosity is represented by 

Ht = Repavt(7 k ) b 

In the inner layer, 0 < y < y c , the expressions for v and £ are 


(*). 


where to is vorticity, and where 


= 0.4y 


y + = v 


:r = «|W| 

Rep\u n 
P 


( 12 ) 

(13) 

(14) 


In the outer region, y > y c , the expressions for v and £ are 

[U m ax ~ Umrn ) 2 


[v) outer = minimum of^J mox , 

[fyouter — VmaxU BL 


4J„ 


(15a) 
(15b) 

where U is the absolute value of velocity, and where the value of Cbl is 1-6. Instead of using 
equation (15a), we have used the following expression 


[v] outer = minimum of 


(' 


{Ut...,, ~ Unun)'- 

3 max 


(15c) 


which has been used by Baldwin since the publication of reference 2. The function J(y) is defined 
to be 


?(y) = yM 


1 — exp 


(-fi)l 


(16) 


and the quantity 7 max is the maximum value of 7(y) that occurs in this equation, and y m ax is 
the value of y at which it occurs. 


Baldwin and Lomax (ref. 2) simulated the effect of intermittency by using 

o — l 

i / r,n \ ' 

Mb = 


i + 5.5(^y 

\ Umax ) 


(17) 


where Ck = 0.3. Further, they simulated the effect of transition by setting pt equal to zero 
everywhere in a normal profile for which the maximum tentatively computed value of pt from the 
foregoing relations is less than a specified value. They have recommended the following condition 

Pt — 0 if [pt) max < -- 14/Zoo 

Visbal and Knight (ref. 19) and York (ref. 20) evaluated the Baldwin-Lomax turbulence 
model for two-dimensional turbulent boundary-layer computations. These investigators made 
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a comparison with the Cebeci model. Based on these comparisons, they recommend that the 
constants Cbl and C k should be, respectively, equal to 1.2 and 0.646 for M <*, = 0.01. This value 
of Cbl would give a lower value of skin friction than that with the value of 1.6 used in this study. 


The Baldwin-Lomax turbulence model is deficient in three aspects, and it gives erroneous 
results for one situation. First, it does not simulate the transition process. There is an abrupt 
change in p-t from a value of zero in laminar flow to a nonzero value in turbulent flow. This 
is corrected by including the transition factor given in equation (8). Second, the model is not 
valid in the far-wake region since the inner eddy-viscosity formulation extends to large values of 
y; therefore, it dominates the outer or wake eddy-viscosity formulation. In the present study, 
this anomaly has not been corrected, since the grid system is not appropriate in the TLNS 
computations. And third, the model does not account for the effect of pressure gradient. This 
deficiency is also not eliminated. The model erroneously gives pt = 0.0 in the wake region when 
the wake is symmetric, because the Van Dreist damping function, equation (14), depends on the 
vorticity on the centerline of the wake, which has the value of zero. This error has been corrected 
by the expression 


(Ht)b = RepavLi tr {~i k ) b 


with the length scale equation now given by 




t nner 


= 0.4y 


with 


_ / 0 . 0 , 

" l (*“ 


1 — exp 


x te)i 


(_y^_ _ x + V 

V 26 26 ) 


if x< 1.0 
otherwise 


(18) 


(19) 


In the outer region, the expressions for v and t are the same as those given in equations (15b) 
and (15c), but the definition of 7{y) is replaced by the following: 


7{y) = y|w 


1 — exp 


(_ _ x^y 

V 26 26 ) 


( 20 ) 


Transition Criteria 


There are a few empirical methods for computing the location at which a laminar flow begins 
to undergo transition to a turbulent flow. The transition location is fixed by experimental obser- 
vation, by laminar separation, or by Michel’s correlation method (ref. 21). Michel’s method is 
not valid for flows with separation; instead it is assumed that the computed laminar separation 
location is the transition location. Michel’s method is embodied in the equation of reference 22, 
that is, 

■ 22 - 400 ' 


1 + 


Re 


ZtT 


Re 


0 46 
Xtr 


( 21 ) 


where subscripts 6 and x designate Reynolds numbers based on momentum thickness and surface 
distance, respectively. 
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Various Models 


The turbulence models and transition criteria have been used in the combinations of table 
1. For example, Model III is a combination of the wake formulation of model (/x t ) a and of the 
boundary-layer formulation of model (fit)b with the empirical transition criteria. Models I- IV are 
used in the IBL method, whereas only Models II and V are used m the TLNS method. It was not 
feasible to use model (nt)a and Model III with the latter method because of large fluctuations in 
the location of the edge of the boundary layer during the iterative solution procedure. 

TABLE 1. - COMBINATION OF TURBULENCE MODELS AND TRANSITION CRITERIA. 



Note: EMP. = empirical; EXP. = experimental; B. L = boundary layer. 


DISCUSSION OF RESULTS 


We have chosen a flow over a NACA 0012 airfoil section for which to compare the TLNS and 
IBL computations with measurements, and with each other. Two entirely different experiments 
at almost the same Reynolds number are chosen in order to compare not only the computations 
with data, but also to compare data of one experiment against the other. The experimental 
data, reported by Gregory and O’Reilly (ref. 23), were obtained during their Session III at 
Re = 2.88 x 10 6 and Moo = 0.16 on a clean (smooth) airfoil, and those reported by Loftin and 
Smith (ref. 24) were determined at Re = 3 x 10 6 and Moo = 0.07. 

The IBL computations are done at Moo = 0.0 assuming that the fluid is incompressible, 
whereas the TLNS computations are carried out at M c 0 = 0.1. The latter computations could 
have been done at the experimental Mach numbers, but this was not done; the reason it was not 
done was the desire to have the same Mach -number influence in both sets of TLNS computations 
for comparison with the IBL computations, and to keep the Mach number low enough for the 
TLNS computations without requiring excessive computational times. Anyway, the computa- 
tional Mach numbers are low enough for compressibilty to have little influence, even at maximum 
lift. This can be verified from the computed Mach-number contours shown in figure 1. Most of 
the compressibility effect is near the nose of the airfoil, where the maximum local Mach number 
is 0.35. 
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Figure 1. - Mach-number contours: NACA 0012, Re = 2.88 x 10 6 , a = 16°, and — 0.1. 


Computational Grids 

. i 

The total number of grid points used in the TLNS calculations is based on the decision to 
use about 85% of in-core memory of one processor of the CRAY X-MP 22 computer. These 
calculations made use of a C-type-grid topology with 257 x 57 node points, and with 36 points in 
the wake region of the airfoil along the free-stream direction. It was generated with an algebraic 
grid generator called CGUAVA based on a procedure formulated by Eiseman (ref. 25). Except 
near the leading edge and the trailing edge of the airfoil, the first grid point off the airfoil surface 
is at a distance of 0.00001L, where L is the chord length. For fully developed, attached turbulent 
boundary layers, this distance corresponds to y + at the first grid point off the surface in the range 
of 0.15 to 3.5, depending on chordwise location and angle of attack, at Re = 2.88 x 10 6 . The first 
grid point upstream of the leading edge and downstream of the trailing edge is at a distance of 
0.001L. The outer boundary is kept at 10 chord lengths from the airfoil. Figure 2 shows a part 
of the grid system. 

The IBL calculations were performed with approximately 80 x-stations on the upper surface 
and approximately 60 on the lower surface, and with 37 y-stations for small angles of attack and 
50-70 for large angles of attack; the y-stations were on the upper surface from the leading-edge 
region to the trailing edge of the airfoil. There are 30 x-stations in the wake, which extends up 
to 3 chord legths downstream of the trailing edge. 
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Figure 2. - C-grid around the NACA 0012 airfoil. 
Transition Locations 


Figure 3 shows the transition locations 1 computational, experimental, and empirical locations 
based on Michel’s correlation method. Here lower surface transition locations at positive inci- 
dence are plotted as if they were upper-surface locations at negative incidence. The transition 
locations are fixed for Re = 2.88 x 10 6 as follows: When no transition bubble is experimentally 
observed, both IBL and TLNS computations have used the experimental transition locations. 
When this bubble is observed, the IBL computations have considered transition location at the 
experimental laminar separation point, whereas the TLNS procedure has taken this location to 
be approximately at the streamwise center of the bubble. The actual computational transition 
locations are slightly different from these locations owing to the fact that there may not be a 
grid point at the experimental locations. In addition to these computations, computations have 
been made with the IBL method using the empirical locations. On the other hand, at a Reynolds 
number of 3.0 x 10 6 , both procedures have used empirically computed transition locations. Again, 
the TLNS computational locations are slightly different owing to the fact that there may not be 
a computational grid point exactly at the empirical location. 

As can be seen from figure 3, for Re = 2.88 x 10 6 the differences in the locations of transition 
arrived at by experiment and the previous procedures were largest, approaching 10% of chord, at 
low angles of attack. On the upper surface, at these angles, the influence of transition is small. The 
discrepancies associated with angles of attack greater than 9° stem from the different approaches 
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Figure 3. - Transition locations on the NACA 0012 airfoil. 


to the transition location in the presence of a separation bubble; the largest discrepancy is less 
than 1% of chord, which is within measurement precision. Furthermore, on the upper surface of 
the airfoil, empirical locations are upstream of the experimental transition locations or laminar 
separation points, except at angles of attack of 4° and 6°. For Re = 3 x 10 6 T the empirical 
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locations on the upper surface of the airfoil are downstream of those for Re = 2.88 x 10°. This 
should not be the case, and at present we have no explanation for it. 


Comparisons of Loads 


The calculated and measured values of lift coefficient Cl are shown in figure 4 as a function 
of angle of attack, and they are in close agreement up to 12°. At a = 16°, the 1BL apprach 
did not work with Model IV, but it did work with Model I (tables 1-3, and fig. 4). The only 
difference between these two models is the specification of the transition location, being, for 
example, at 0.8% of chord for Model IV and at 0.03% for Model I at Re = 2.88 x 10 6 . At 
higher angles of attack, the expected fall off in Cl does not occur. However, as it is shown in 
reference 1, a small change in the location of transition can cause lower values of calculated Cl- 
The magnitude of the required change in the location of transition is similar to that which can 
be resolved experimentally. The TLNS approach cannot predict this effect, probably because of 
the following reasons: (l) it does not use a time-accurate numerical scheme, a requirement for 
predicting stalling characteristic, which is by definition unsteady; and (2) the normal pressure 
gradient at the surface of the airfoil is assumed to be negligible. On the other hand, the IBL 
approach does a better job of predicting this lift fall off than the TLNS approach. This conclusion 
is misleading because, if the wake curvature effect were included in the IBL approach then the 
lift coefficient would be slightly higher than that computed. 

TABLE 2. - LIFT AND DRAG COEFFICIENTS 
FOR THE NACA 0012 AIRFOIL AT Re = 2.88 x 10 6 . 


a, 

deg 


C 

L 



c 

D 


EXP. 

IBL 

TLNS 

EXP. 

IBL 

TLNS 

1 

IV 

V 

1 

IV 

V 

0 


00 

0.0 

0.0 


0.00596 

0.00562 

0.00568 

2 


KEH 


mam 


0.00611 

0.00571 

0.00572 

4 

0.440 

0.421 

0.418 

Ell 


0.00633 

0.00689 

0.00750 

6 

0.650 

0.639 

0.637 

Ell 

1 

0.00689 

0.00781 

0.00921 

8 

0.850 

0.873 

0.856 


0.0112 

0.00899 

0.00921 

0.01207 

10 

1.055 


1.081 

mm 

0.0134 

0.01070 

0.01136 

0.01547 

12 

1 240 

1.256 

1.268 

E9 


0.01294 

0.01306 

0.01962 

14 

1.400 

1.422 

1.448 


0.0244 

0.01662 

0.01447 

0.02469 

16 

1.490 

1.545 


19 

0.0338 

0.02315 


0.03121 


The influence of the turbulence model was also investigated, and the results are presented in 
tables 2 and 3. The effect on Cl for Re = 2.88 x 10 6 is negligible, except at a = 16°, as discussed 
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Figure 4. - Variation of Cl, with a for the NACA 0012 airfoil. 


above. The lift coefficients for Re = 3.0 X 10 6 are nearly the same up to a = 10° for turbulence 
Models I, II, and III. When the angles of attack are greater than 10°, the IBL approach does not 
converge with Model II, and lift values obtained with Model III are not satisfactory The TLNS 
approach does not converge with Models I and III primarily owing to the numerical scheme that 
is not time-accurate. 

TABLE 3. - LIFT COEFFICIENTS FOR THE NACA 0012 AIRFOIL AT Re = 3 0 x 10 6 . 


a, 

deg 

c L 

EXP. 

INTERACTIVE 

TLNS 

1 

II 

III 

II 

0 

0.03 

0.0 

00 

00 

0.0 

2 

0.23 

0.2097 

0.2130 

0.2131 

0 2166 

4 

0.44 

0 4202 

0.4265 

0 4266 

0 4310 

6 

0.64 

0 6372 

0 6472 

0 6488 

0 6450 

8 

0 88 

0 8726 

0.8832 

0 8818 

0.8599 

10 

1 10 

1 0775 

1.1110 

1.1102 

1.0614 

12 

1.26 

1 2665 


1 3143 

1 2608 

14 

1 43 

1.4293 


1 5190 

1.4489 

16 

1.52 

1.5449 


1.7163 

1.6313 


As shown in figure 5, a comparison of the NACA data (ref. 24), the NPL data (ref. 23), and 
the NASA data (ref. 26) in terms of the lift-curve slope using Kaplan’s rule (ref. 27), raises some 
question concerning the NACA data. Both the IBL and TLNS procedures give lift-curve slopes 
that compare well with the NPL and NASA data. 

The drag results of tables 2 and 4, and figure 6 require more detailed discussion. The IBL 
approach has used the velocity defect in the wake to compute the drag coefficient Cp, whereas, 
the TLNS procedure has integrated forces around the airfoil surface. The pressure drag ranges 
from about 13% to 86% of the total drag, for angles of attack of 0° to 16°. It is likely that 
the pressure drag integrated using surface-pressure distribution is determined with acceptable 
accuracy, especially in view of the quality of the calculations suggested by the comparisons of 
computed and experimental pressure coefficients discussed below. At zero angle of attack, the 
total drag computed by the IBL method is within two counts in hundred thousandths and in ten 
thousandths, respectively, for Re = 2.88 x 10 6 and Re = 3.0 x 10 6 , of that computed by the TLNS 
method; and both values agree well with the experimental values. At angles of attack of 0° and 2°, 
the total drag computed by the IBL method is within six counts (hundred thousandths) of that 
computed by the TLNS method, with both methods utilizing experimental transition locations; 
values of both methods agree well with the experimental values. At higher angles of attack, before 
the appearance of the transitional bubble, the computed results of the IBL method are lower than 
those determined by experiments and the TLNS method. At still higher angles of attack in the 
presence of this bubble, with turbulent flow separation and lower pressure coefficients along the 
aft region of the upper surface, the IBL values are still lower. This drop is probably a result of the 
fact that the wake-curvature effect and the cross-stream pressure gradient effect are not accounted 
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Re X 10* 6 

Figure 5. - Lift-curve slope for the NACA 0012 airfoil. 

for in this approach. The cross-stream pressure gradient increases the momentum thickness of 
the wake, and consequently, the total drag. This conjecture is consistent with the findings of 
reference 6, in which the use of the present ( fi t ) a formulation for wake flows was examined. 

The influence of the turbulence model as indicated in tables 2 and 4, although more pronounced 
in Cp than in Cl, is unlikely to explain the discrepancies of figure 6. If the values of constants 
Cbl and Ck recomended in references 19 and 20 were used in the formulation of [fit) a, then the 
friction drag coefficients would be lower than those computed with the present formulation. 

I 

Further examination of the experimental values of Cd showed that they are different from 
those reported by Loftin and Smith (ref. 24) for Re — 3.0 x 10 s . For example, at an angle of attack 
of 12° the latter NPL reference has reported Cd higher by 14%. Since the NPL data are dated 
25 years later than those of reference 24, they are probably more accurate. The TLNS procedure 
appears to do a better job of predicting drag than the IBL procedure, which underpredicts at 
high angles of attack. However, it is evident that the uncertainties associated with the calculation 
procedures and with measurement of Cd require further examination. 

The experimental moment-characteristic documented in the NACA report is quite different 
from that documented in the NPL report (fig. 7). Assuming that the latter data are accurate, 
the TLNS procedure does not well predict the moment coefficient at high angles of attack. 

A study was conducted to provide some indication of the influence of the outer boundary 
conditions on loads acting on the airfoil using the TLNS approach. Two different locations are 
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TABLE 4. - DRAG COEFFICIENTS FOR THE NACA 0012 AIRFOIL AT Re = 3.0 x 10 6 . 



C D 

a, 

deg 

EXP. 

INTERACTIVE 

TLNS 

1 

II 

III 

II 

0 

0.0059 

0.00585 

0.00592 

0.00595 

0.00614 

2 

0.0062 

0.00606 

0.00608 

0.00619 

0 00632 

4 

0.0070 

0.00635 

0 00645 

0 00648 

0 00734 

6 

0.0084 

0.00700 

0.00704 

0.00700 

0.00885 

8 

0.0100 

0.00881 

0.00892 

0.00899 

0 01227 

10 

0.0130 

0 01062 

0.01065 

0.01071 

0.01585 

12 

0.0158 

0.01260 


0.01261 

0.02005 

14 

(0.0204) 

0.01638 


0.01495 

0 02567 

16 

(0.0264) 

0.02315 


0.01853 

0.03206 


Note: The values in parenthesis are determined by extrapolation. 


considered for the outer boundary: 10 chord lengths and 20 chord lengths from the the trailing 
edge. The increase in the distance from 10 to 20 chord lengths is achieved by adding grid points, 
resulting in a 273 x 65 grid-system with 44 grid points in the wake, such that the grid distribution 
in the vicinity of the airfoil is essentially identical for these two different locations of the outer 
boundary. Further, computations are carried out with and without the circulation correction (ref. 
11) to the free-stream values, as discussed in the appendix. The results of this study are shown in 
figure 8 for a = 2°. The variations in the lift, drag, and moment coefficient are, respectively, 1%, 
6 counts (hundred thousandths) or 1%, and 4.5%. These variations for lift and drag coefficients 
are within the experimental accuracy. Therefore, all results presented in this report, except those 
mentioned in figure 8, have been obtained with a 257 x 57 grid-system, with the outer boundary 
at 10 chord lengths, and without the circulation correction 


Some Flow-Field Details 


Distributions of pressure coefficient are shown on figures 9-12. Figures 9 and 10 compare the 
TLNS and IBL results, respectively, with the experimental data of reference 23 for a range of 
angles of attack; and figure 11 compares the results of the two procedures for angles of attack 
of 10° and 12°. As can be seen, the calculated results of the two procedures agree very well 
with each other and also with the experimental data. In figure 9, at 14° and 16°, the TLNS 
results show some irregularities near the suction peak. This is caused by improper transition 
location. This conclusion is supported by the results shown in figure 12. There is no irregularity 
at the above angles of attack, but there is at the angle of attack of 6°. The latter is caused by 
the computational transition location being downstream of the possible experimental location. 
Such an irregularity is not there at the corresponding angle of attack at Re = 2.88 x 10 6 (fig. 
9). Observe that the empirical transition location at Re = 3 x 10 6 is downstream of that at 
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Figure 6.- Drag characteristic of the NACA 0012 airfoil. 



MOMENT COEFFICIENT, C M 


.02 r 




(a) Re = 2.88 x 10 6 (b) Re = 3.0 x 10 6 


Figure 7. - Moment coefficients for the NACA 0012 airfoil. 


Re = 2.88 x 10 6 , both being downstream of the experimental location at the lower Reynolds 
number. 


O A □ WITHOUT CIRCULATION 



RECIPROCAL OF DISTANCE, 1/R 

Figure 8. - Effect of outer boundary conditions on loads: 

NACA 0012, Re = 2.88 x 10 6 , = 0^ a = 2°. 

04 

Pressure contours around the leading edge and the trailing edge of the airfoil at a = 2° and 14° 
are shown in figure 13. The normal pressure gradient near the leading edge is negligible at small 
angles of attack, but not at large angles of attack. For the range of angles of attack investigated, 
immediately downstream of the trailing edge the normal pressure gradient is negligible. These 
conclusions arc confirmed by the streamline plots. As shown in figure 14, the streamline curvature 
effects are expected to be important at high incidences near the leading edge, and in the wake 
region where the flow tends to go back to the free-stream direction. 

Calculated distributions of displacement thickness 6 ’ are reproduced in figure 15. Corre- 
sponding experimental data are not available, but the results are interesting in that they show 
the rapid growth of 6* with angle of attack; and both numerical schemes give similar results. 
Because of the uncertainity of defining the edge of the boundary layer with the grid system used 
near the aft portion of the airfoil in the TLNS method, the smaller of (l) the inviscid velocity at 
the surface of the airfoil, based on the Bernoulli’s equation and the ideal gas relation, and (2) the 
maximum velocity along the grid line emanating from the surface, is taken as the edge velocity 
for computing 6 ’. Consequently, there is a kink in the 6’ distribution when a switch is made 
from one family of edge velocities to the other. 
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Figure 9. - Comparison of pressure distributions: TLNS results and experimental data 

at Re = 2.88 x 10 6 . 


The TLNS method is able to compute transitional separation bubbles. With a Navier-Stokes 
procedure, this is the first time such a bubble has been predicted as a part of the overall prob- 
lem of computing the flow past an airfoil. Figure 16 shows that the computed separation and 
reattachment points do not agree with the experimental points. In addition to this, the compu- 
tation shows a transitional bubble at a — 8°, whereas the experiment does not. These differences 
are owing to the lack of grid resolution, to imprecise computational transition location, and to 
experimental inaccuracies associated with a small transition bubble. 

The influence of transition location on a transitional separation bubble is illustrated in figure 
17. If the transition location is upstream of the location where it should be, the transition bubble 
cannot be observed in computation (curve (a) of fig. 17). On the other hand, if it is downstream of 
the location where it should be, it is difficult to get a converged solution, supporting a conclusion 
by Cebeci et al. (ref. 1) and by Carr and Cebeci (ref. 31) from calculations based on the IBL 
method. Cases (c) and (d) in figure 17 are solutions at two different iteration numbers with the 
transition location being the same. It was not possible to obtain a converged solution. Cases (a) 
and (b) are converged solutions. 

Turbulent separation regions from the TLNS method are reported in figure 18. Note that 
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Figure 10. - Comparison of pressure distributions: IBL results and experimental data 

at Re = 2.88 x 10 6 . 

in this figure the vertical scale is 20 times larger than the horizontal scale. Experimental and 
computed separtion points from TLNS and IBL computations do not agree, probably owing to a 
lack of grid resolution, turbulence model, and experimental inaccuarcies. 


Computational Efficiency 


The evaluation of the two procedures can be made, in part, by information of the type referred 
to in the previous paragraphs. It is often necessary also to consider the cost and ease of achieving 
the results. In this respect, it is important to note that one iteration of the IBL procedure 
and that of the TLNS procedure require about 0.1 min of CPU time on an IBM 3381 computer 
and about 0.005 min on one processor of a CRAY X-MP 22 computer, respectively. The latter 
procedure, which is not developed for computing flows of essentially incompressible fluid , is quite 
costly. For the IBL method, the convergence criterion is 0.1% residual in the conservation of 
momentum, requiring up to approximately 10 iterations; whereas, for the TLNS methods, the 
convergence criterion is 10“ 11 or lower residuals in the governing equations, requiring several 
thousand iterations. 
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Figure 11. - Comparison of pressure distributions: TLNS and IBL results. 


Different criteria are available for terminating an iterative procedure. They are based on 
either computational cost or computational accuracy. The primary consideration in the TLNS 
computations has been accuracy. Criteria based on accuracy vary from plotting or engineering 
accuarcy to computational machine accuracy, with the former being subjective and the latter being 
machine dependable. This is illustrated with the convergence histories for the TLNS computation 
shown in figures 19 and 20. 

The convergence histories of Cl and Cp shown in figure 19(a) suggests that the converged 
solution is achieved within 1,200 iterations. But, the machine accuarcy is not achieved even after 
31,000 iterations (fig. 19(d)). One percent error in Ci is achieved long before a corresponding 
error in Cp. An error of less than one count (ten thousandths) in Cp is reached in 14,000 
iterations. On the other hand, the convergence history of the L 2 -norm of the residual of the 
governing equations, taken together or individually, would lead to somewhat different conclusions. 
An example is presented in figure 20 for the mass conservation equation. When this equation is 
satisfied its residual is zero. If the convergence criterion is that the residual should decrease by 
four orders of magnitude, then in about 1,200 iterations a converged solution is obtained, whereas 
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CHORDWISE DISTANCE, x + 0 05 a° I 

Figure 12. - Pressure distributions determined with the TLNS method and Model II 

at Re = 3.0 x 10 6 . 

the convergence criterion of six orders of magnitude results after about 6,500 iterations (fig. 
20(a)). A further 15,000 iterations hardly not reduce this residual, with the residual fluctuating 
around 10 -1 ° for the first 7,500 iterations (fig. 20(b)) and around 2 x 10 -u for the second 7,500 
iterations (fig. 20(c)). Eventually, the residual does decrease monotonically as seen in figure 20(d). 
Again, more than 31,000 iterations are required for machine accuracy. In these computations, 
the reference time-step (eq. 45) was varied from 1 to 10, as it was not possible to maintain the 
same value for avoiding divergence. 

The large number of iterations required for the TLNS is a result of the following: (1) the 
numerical scheme is not linearized with respect to the effective viscosity, and (2) the numerical 
scheme is not efficient for low Mach numbers. The latter is substantiated by the convergence 
histories shown in figures 21 and 22 for M 0 0 = 0.5. A possible modification of the present scheme 
for efficient computations for low Mach number is discussed by Turkel (ref. 32). 


UPSTREAM INFLUENCES 


The boundary-layer theory is based on the assumption that the shear layer is thin and on the 
condition that this layer grows only slowly in the general direction of flow. As a consequence of the 
above assumption and condition, the upstream influence within the thin layer is not accounted 
for by the first-order theory and is only partially accounted for by the second-order theory. This 
is the fundamental limitation of the boundary-layer theory. | 
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Figure 13. - Pressure contours around the nose and the tail of the NACA 0012 airfoil: 

Re = 2.88 x 10 6 and Moo = 0.1. 
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Figure 16. - Transitional-separation bubbles determined with the TLNS method and Model V 
on the NACA 0012 airfoil: Re — 2.88 x 10 6 and M <*> = 0.1. 
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Figure 17. - Effect of transition location on transitional-separation bubbles: 

Re = 2.88 x 10 6 , = 0.1, a = 10°. 

| 

With s as the streamwise coordinate, a shear layer, whose thickness is 6 at a distance L from 
its origin, is considered to be thin if 36/ ds is of order 6/L and is much smaller than unity, that 
is, if 36 /ds ^ ^ 1. Note that it is the local rate of change that matters. In this case, the ^Javier— 
Stokes equations reduce to the boundary-layer equations. These latter equations are useful only 
if they are in a coordinate system such that s is in the direction of the thin shear layer. When this 
thin-shcar-laycr approximation is applicable, the ratio of neglected to retained stress gradients is 
0\{36 / ds) 2 \ for laminar flows and 0(36/3s) for turbulent flows. Since ( 36/ds ) is generally larger 
in turbulent flows than in laminar flows, this approximation is less accurate in turbulent flows. 
Typically, the neglected streamwise viscous diffusion and turbulent diffusion are, respectively, 
0(5.3Re ~ l ) and Oft^Re- 1 / 5 ) for Re < 10 7 . 

The upstream influence is provided by the pressure field in subsonic regions, by transport of 
momentum in the upstream direction in separated flows, and by streamwise gradients of viscous 
and turbulence stresses, as indicated by the direction of the arrows in figure 23. The first or 
the last of these influences makes the flow field and the governing equations elliptic in the spatial 
dimensions. The correction owing to the displacement effect can indirectly transmit the upstream 
influence of the pressure field via the external flow, but not through the shear layer. A correction 
owing to displacement thickness in the “inviscid” flow at a station leads to a correction in the 
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Figure 18. - Turbulent separated regions obtained with the TLNS method and Model V 
on the NACA 0012 airfoil: Re — 2.88 x 10 6 and Moo = 0.1. 


pressure distribution upstream of this station. This procedure is a part of the IBL theory. On 
the other hand, the correction owing to curvature effect accounts for the variation of the pressure 
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Figure 22. - Convergence history for mass conservation equation: 
Re = 6.0 x 10 6 , a = 2°, and = 0.5. 




field across the shear layer, but not for the upstream influence of this field through this layer. In 
addition, this correction is limited to shear layers that are thin, since the curvature of streamline 
varies substantially across thick shear layers. Usually, in high-i?e, laminar separated flows there is 
significant upstream influence of the velocity field relative to that of the pressure field within the 
shear layer so that the latter can be neglected. The boundary-layer theory does not account for 
the upstream influence of streamwise gradients of viscous and turbulence stresses. These stresses 
are generally negligible at normal (high) Reynolds numbers. 


INVISCID 

BOUNDARY 



The main difference between the present IBL procedure and the TLNS procedure is iden- 
tified by considering what upstream influences these procedures account for. Neither handles 
the upstream influence of viscous and turbulent stresses. Both handle the upstream transfer of 
momentum. The difference between these procedures arises when the upstream influence of the 
pressure field is dealt with. In order to explain this, the upstream influence of the pressure field 
is broken up into two parts: one through the shear layer and the other through the inviscid re- 
gion. The TLNS formulation accounts for both of these parts; but the IBL formulation accounts 
only for the second part, through the consideration of the displacement thickness. This is the 
fundamental difference between these two formulations. 

When a shear layer is neither thin nor slender, the Navier-stokes equations are necessary. The 
upstream influence of both the pressure field and streamwise gradients of viscous and turbulent 
stresses are then important. The upstream influence of the pressure field is much more pronounced 
than in a thin shear layer, owing to significant normal pressure gradients; for example, in large 
separated regions, in the trailing-edge region of an airfoil, and in interactions of a shock wave 
and a boundary layer. Further, the streamwise gradients of stresses cannot be neglected, for 
example, in low-Reynolds-number flows and when the velocity field encounters rapid variations in 
streamwise direction, such as flow past the trailing edge of an airfoil, across a shock wave within a 
boundary layer, and around a point of flow separation. It remains to establish the magnitude and 
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significance of error introduced by the IBL approach vis-a-vis the TLNS approach in the above 
situations that are of relevance in high-speed aerodynamics. 


CONCLUDING REMARKS 


Numerical procedures for the solution of the TLNS equations and for the interaction of solu- 
tions to inviscid and boundary-layer equations were described and evaluated. To allow appraisal 
of the numerical and fluid dynamic abilities of the two schemes, they were applied to one air- 
foil over a range of angles of attack. Calculations were performed with algebraic eddy-viscosity 
formulations, and consideration of transition was included. In general, the IBL procedure is 
computationally efficient but theoretically restrictive relative to the TLNS procedure. 

Some specific conclusions can easily be extracted from the preceding section. It is clear, for 
example, that both procedures will allow the calculation of lift coefficient with acceptable accuracy 
at angles of attack that cover laminar separation bubbles but that fall short of Ci max . The TLNS 
procedure appears to better predict drag than the IBL procedure. The latter procedure can be 
improved by including the cross-stream pressure gradient, and the wake streamline curvature 
effect to better predict the drag force. 

Based on the TLNS results, the far-field wake turbulence model seems to have negligible 
effect on the loads acting on the airfoil. In common with the findings of references 6 and 34, the 
calculations of pressure and lift confirm that the algebraic eddy-viscosity formulations adequately 
represent the boundary-layer flows and the wake up to angles of attack close to that of Cx, max , 
as long as the turbulent separated region is thin. 

The location of transition is crucial. Either some significant feature may be missing in a 
computation or an iterative procedure may not converge, depending on whether the location in 
the computation is upstream or downstream of that feature in reality. Further work is necessary 
to determine the role of transition at high angles of attack, and to confirm and extend this report 
and the reports of Cebeci et al. (ref. 1) and of Carr and Cebeci (ref. 31) concerning the influence 
of the transition location on the flow properties and the break down of calculation procedure. 

The studies in reference 1 indicate that the shape of the airfoil also plays an important role 
in calculating its lift and drag coefficients. Additional calculations, similar to those discussed in 
this paper for a symmetrical airfoil, should be performed for cambered and supercritical airfoils. 

When the upstream influence of pressure field through the shear layer is important, the TLNS 
procedure has an edge over the IBL. There is no doubt that the TLNS equations are more general 
than those associated with the IBL procedure, even if the IBL equations are modified to include 
the cross-stream pressure gradient. On the other hand, the cost of using the computer program 
that embodies the TLNS equations is considerably greater than that of the program that embodies 
the IBL equations. This cost may limit the use of the TLNS equations in a design cycle with 
presently available computing power even for two-dimensional flows. Therefore, it is necessary 
to improve the numerical algorithm of the Navier-Stokes procedure and, in the meantime, to 
attempt to establish those flows for which its added generality is required. 
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APPENDIX 


THE TLNS METHOD 


I 


j 

I 


The derivation of the approximate-factorization scheme presented below is different from that 
given by Beam and Warming (ref. 3) for the Navier-Stokes equations in Cartesian coordinates and 
extended by Steger (ref. 4) for these equations in curvilinear coordinates. The steps involved in 
formulating the present scheme are as follows. First, the governing equation is evaluated at a new 
or unknown time-level (n + 1) and the time-derivative is approximated by a backward difference 
formula. Second, assuming that the effective viscosity and thermal conductivity are not functions 
of time, the spatial derivatives at the (n + 1) level are expressed by a Taylor series expansion 
in terms of known quantities at the ( n ) level, and of time-derivatives of dependent variables at 
the (n + 1) level. This leads to linearization of the nonlinear terms in governing equations. The 
Jacobian and metric coefficients of spatial coordinates and the time- variable are considered to be 
known at the (n+1) level. Third, the left-hand side of the equation resulting from the previous step 
is approximately factored to form an ADI procedure. Fourth, the continuous space derivatives are 
replaced by finite differences. Fifth, numerical dissipation is added in the numerical procedure. 
Sixth, physical and analytical boundary conditions are translated into finite-difference boundary 
conditions. Instead of formulating these boundary conditions implicitly and incorporating them in 
the ADI procedure, boundary values of the dependent variables are determined using the previous 
time-step solution. Seventh, the solution of the equation formed in the fifth step determines the 
time-derivative of the dependent variable at the (n+l) level. This time-derivative then is replaced 
by a backward finite-difference formula in order to determine the dependent variable at the (n + l) 
level. 

The above steps are explained below in detail for the thin-layer Navier-Stokes equations. Note 
that the derivation of the approximate factorization scheme for the Navier-Stokes equations follow 
the same steps. 


Linearization and Approximate Factorization 
The nondimensional form of the TLNS equations are written in the following form: 

dQm 


dr 


= 9?r 


with 


where 
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1 dT; 


1=1 


Re di 2 


X 


= E (w ^ 

i=i ' 


dxj 


( 22 ) 


38 



The Cartesian, conservative dependent variables, convective and pressure terms, and diffusive 
terms are associated with Q m , C mt , and respectively. The Jacobian of transformation from 
the Cartesian coordinates x x to curvilinear coordinates £, is denoted by V. The subscript m is 
the equation number, varying from 1 to 4. 

Let the numerical solution of equation (22) with appropriate boundary and initial conditions 
be known at time- level (n). The numerical solution of this equation is required at time-level (n+l) 
by using an implicit scheme. Therefore, this equation is evaluated at the ( n + 1) time-level, and 
the time-derivative is replaced by a backward difference formula, 




+ 0(Ar p ) 


(23) 


where 6 is a backward finite-difference operator that is not of the Pade type. This backward 
difference formula is not presented at this point because it is not required in the next few steps. 
It is defined in the seventh step below along with the value of p. 

After substituting formula (23) for the time-derivative in equation (22), we have 



»+ 1) + 0(Ar p ) 


(24) 


We have unknown quantities on both side of equation (24), making this equation and the numerical 
algorithm developed below implicit. 

The terms represented by 9?^ +1 ^ on the right-hand side of equation (24) are expanded in a 
Taylor series of time r. The expansion of C mt -term is 
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where 
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(26) 


and equation (23) are used. The vector of the viscous term contains first derivatives. These 
are derivatives of functions that are of the form (3 = (3(Q ). Symbolically, this is written as 
= i,/? 2 ,...) = T^,(/?). Instead of expanding directly (ref. 3) in a Taylor series, it is 

indirectly expanded by writing each /? in a Taylor series: 


P n +i = 
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The viscous term at the (n+l) level is the expressed as 

MM 


^ ^m(/?) y w+1) 


ae 2 


1 ( * + Ar a ^^ (n + ! ) ) 1 + 0( A r 2 ) + 0[ArP +1 ] 




(28) 


39 



where 
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By using equations (25) and (28), we can reformulate equation (24) to 
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= 9?^^ + 0(Ar 2 ) + 0[At p ] 

where the Jacobian and metric coefficients are evaluated at the ( n + 1) time-level in 9?™*^. This 
is the only difference between this latter quantity and The approximate factoring of the 

left-hand side of the above equation and the omission of the truncation error on the right-hand 
side lead to the following alternating direction algorithm in “derivative” form: 
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(31) 
in the fi- 


with a factorization error 0(Ar 2 ). The first factor is in the ^-direction rather thanj 
direction so that the effect of the boundary conditions is transmitted quickly from the surface 
of the airfoil to the interior of the flow field. The Roman numerals are used to indicate dummy 
time-derivatives . 


Numerical Dissipation 

Minimum numerical dissipation is used so that stable numerical solutions are obtained in 
reasonable amount of computer time. The added dissipation is patterned after that developed by 
Jameson et al. (ref. 12). The form of this dissipation and the various constants appearing in it 
are kept the same for all cases. The numerical dissipation D m is given by 
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where <p a = . The subscript a is either 1 or 2. In the present study, ei and e- 2 are, respectively, 

0.64 and 0.07. The term is defined as 
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where U and a are the contravariable velocity and the speed of sound, respectively, and where 
the constant c — 64. Equation (33) was specifically formulated for the present study. Observe 
that the numerical dissipation in the ^-direction, the direction in which viscous and turbulent 
stresses are dominant, is much less than that in the other direction. 

The above fourth-order dissipative term damps the short wavelength (2A£ t ) oscillations. Since 
it is of higher order than the second-order-accurate central differences, it does not disrupt the 
formal accuracy of the numerical scheme. Although numerical dissipation is not, in conventional 
terminology, part of the turbulence model, it is in practice. Consequently, the analytical form of 
the turbulence model is not sufficient to describe its effect on an actual calculation. The numerical 
turbulence model includes the internal logic controlling the local evaluation of parameters such 
as mixing length, the discretization error, and the numerical dissipation (ref. 8). It seems that 
the only way the effect of numerics can be separated from the analytical turbulence model is by 
grid-refinement studies leading to elimination of numerical dissipation. In the absence of this, the 
final judgment of the computed results is generally based on a comparison with some experiment 
or theory. 

The numerical dissipation given in equation (32) is treated the same as the right-hand side of 
equation (22). This leads to linearization of D m resulting into two terms, one on either side of 
equation (29). Without the truncation error, this equation becomes 
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where 
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The variable is assumed to be time-independent because its linearization produces terms that 
are difficult to handle and time consuming. The approximately factored form of equation (34) 
leads to pentadiagonal-block systems. 


Spatial Differencing 


Except for the viscous terms, the spatial derivatives appearing in equation (34) are approxi- 
mated by three-point central-difference approximations in the interior of the computational do- 
main. At the boundaries, the metric coefficients are evaluated using two-point sided-differences. 
Vorticity and effective viscosity are computed at the grid-point location of the mesh system. 

The evaluation of the viscous terms requires a special treatment because the effective viscosity 
is a function of spatial coordinates. The form of a typical viscous term and, for example, its 
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second-order finite-difference form (ref. 4), are as follows: 
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Here, the subscript indices denote grid-point locations either in the £j- or £ 2 -direction. 


(35) 


The fourth derivatives in the dissipative terms are evaluated as follows 
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where the various operators are defined as 

^a(;) = {o Sah + l] D b + 1) +0 Uj) D^)) 

= Q m(j ) ~ Qm(j-l) 

A$ u Qm(j) = Qm(j + 1) ~ Qm(j) 


Boundary Conditions 


Boundary conditions for the equations of motion are determined by mathematical and phys- 
ical considerations (ref. 8). The mathematical character of these governing equations dictates 
the number and type of boundary conditions that determine, along with initial conditions, the 
well-posedness of a problem, that is, the existence and uniqueness of solution. Further, this math- 
ematical character is determined by the theory of characteristics. The mathematical character 
of the system represented by the linearized form of equation (22) is incompletely hyperbolic or 
hyperbolic-parabolic. It is parabolic in the (£2 — r ) plane. 

The rigid-wall boundary provides velocity and temperature conditions. The behavior of a real 
gas at ordinary conditions (Knudsen numbers less than 10~ 2 ) is accurately described by the no- 
slip and no-temperature-jump conditions; the latter condition, corresponding to zero temperature 
gradient at the wall, is used. These are the only two physical conditions available. The adiabatic 
condition for internal energy at a solid wall is 



where 1 is internal energy per unit mass, 
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The determination of surface density is carried out with a new governing equation formulated 
by combining the momentum equations to form the normal derivative of pressure, 


dp 

dn 



+ contribution of viscous 


terms 


(38) 


where n represents the direction normal to £ 2 -surface. Since equation (38) is derived from the 
momentum equations, it is not a boundary condition. In the present study, the contribution of 
the viscous terms is neglected and dpj dn = 0. 

At the inflow boundary and at tangent-flow open boundaries, the fluid is assumed to be in- 
viscid. This leads to specification of boundary conditions that correspond to those for the Euler 
equations. These conditions are based on local one-dimensional characteristics. The local char- 
acteristics or eigenvalues determine the number and the admissible forms of boundary conditions 
based on one-dimensional analysis. For a hyperbolic system, the eigenvalues are real. The number 
of negative eigenvalues with distinct eigenvectors determines the number of boundary conditions. 
This number is the same as the number of inward characteristics into the computational do- 
main. In addition, characteristic relations along a characteristic can be extrapolated along this 
characteristic. 


Locally one-dimesional Riemann invariants are given in terms of the normal velocity compo- 
nents as 
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where 

(d£ 2 /dx i) {dti/dx 2 ) 

w n = -u -—V 

[(dti/dx i ) 2 + (d£ 2 /dx 2 ) 2 ] [(df 2 /dxi ) 2 + (3f 2 /3x 2 ) 2 ] 7 

The subscript n designates the direction away from the boundary. The variables u and v represent 
the Cartesian velocity components (nondimensionalized with respect to the free stream velocity). 
These equations along with conditions on the tangential velocity wt and entropy S = ln(p/ p) 1 
determine Q's. The Riemann invariants are associated with the characteristic velocities 


Ai = w n — a and A 2 = w n + a 


When the inflow is subsonic, iv n < 0, and Ai < 0, the variables I\, w t and S are all set to free 
stream values. 


w _ (df 2 /dx 2 ) ^ (gWgfi) v 

[(*fr/3xi) 2 + (d6/dx 2 ) 2 ] 1/2U °° [(dfr/dx!)* + (dt 2 /dx 2 pi 1/2V °° (40) 

S = 1 

where 7 is the ratio of specific heats. When the other characteristic velocity X 2 is greater than 0, 
1 2 is extrapolated from the interior. 
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At the downstream (outflow) boundary, the flow is subsonic and the wake region of the airfoil 
is turbulent. Therefore, fluid density and momentums are extrapolated from the interior and 
energy is computed utilizing the value of free-stream pressure. The wake cut in a C-grid system 
is treated explicitly by specifying the averaged values from either side of this cut. 


When it is stated, the free stream velocities and the speed of sound have been corrected for 
the circulation induced by the airfoil following Salas et al (ref. 11). The corrected free-stream 
velocities are 


«coo = cos(a) + C sin(0) 
v coo = sin(a) - C cos(0) 


(41) 


where 

c 

4nr\l — sin 2 (0 — a)] 

with Ci and a being the lift coefficient and the angle of attack, respectively. A point on the 
outer boundary is fixed by the polar coordinates (r, 0) with the quarter chord location as the 
origin of the coordinate system. In addition, a corrected speed of sound that enforces a constant 
free-stream enthalpy is determined from 


where 
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Numerical Solution Procedure 


(42) 


Because we are interested in only steady-state computations, the numerical scheme given 
by equation (34) is modified in two ways. First, a diagonal form of this scheme developed for 
the Euler equations by Pulliam and Chaussee (ref. 10) is used. (Note that this diagonal form 
decreases the time accuracy of the numerical scheme.) Second, a spacially varying time-step is 
used. Both of these procedures reduce the number of iterations required for converging to a 
steady-state solution. 


The diagonal algorithm reduces the (4 x 4) block pentadiagonal inversion required for solving 
the approximately factored form of equation (34) to (4 x 4) matrix mutiplies and scalar pentadi- 
agonal inversions. This algorithm cuts the operation count of the former procedure by one half, 
and makes it easier to vectorize the computer code. Consequently, the total computational work 
is decreased (refs. 10, 35). 


Because the diagonal algorithm is developed only for the Euler equations, the viscous term on 
the left-hand side of equation (34) is negelected. The flux Jacobian matrices for the convective 
terms in this equation are diagonalized using real eignvalues and a complete set of eignvectors of 
each of these Jacobians (ref. 36): 
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where the eignvectors of the Jacobian matrix are the columns of the matrix r mt . Although the 
latter matrix is a function of spatial coordinates £,, it is assumed not to be. With the above 
assumptions, the approximate factored, diagonalized algorithm of equation (34) is the following 
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(44) 


where N m = T m \T m i. 

The above noniterative, solution algorithm consists of five one-dimensional sweeps. Each 
sweep requires the solution of a linear system involving either a matrix multiple or a scalar 
pentadiagonal system. Except for the last sweep, each sweep determines the dummy value of the 
time-derivative of Q m at the (n + 1) time-level. The last sweep determines its true value. 

The following variable time-step is used to accelerate convergence: 


At = 



(45) 


where the time-step is defined as a local function of the mesh spacing and flow-field variation. 
Thus, an attempt is made to use a uniform Courant number throughout the computational 
domain. 

The finite-difference time-derivative is replaced by a finite-difference formula; for example, the 
following first-order formula is used in the present study: 



The solution of equation (46) determines Q m : 
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To accelerate convergence to a steady state, a good initial guess on a fine grid system, say 
257 x 57, is first obtained. This initial guess is determined by solving the governing equations 
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on a 65 x 35 grid for a few hundred iterations, interpolating the solution, including the eddy 
viscosity, up to a 129 x 45 grid, again solving the governing equations for a few hundred iterations 
on this grid, and then again interpolating the solution and eddy viscosity onto the fine grid. The 
converged solution is then determined on this grid. 
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